home *** CD-ROM | disk | FTP | other *** search
/ Mac-Source 1994 July / Mac-Source_July_1994.iso / C and C++ / Science⁄Math / fastfourierassembly.c
Encoding:
Internet Message Format  |  1994-03-25  |  7.5 KB  |  [TEXT/R*ch]
  1. From: Tom Barrett <barrett@pacific.mps.ohio-state.edu>
  2. Date: Wed, 18 Aug 93 09:35:58 -0400
  3. To: macgifts@mac.archive.umich.edu
  4.  
  5. This is an updated version of /info-mac/sci/fft-in-asm-src.txt
  6.  
  7. * now works for data sets > 64k points
  8.  
  9. * other minor changes for speed improvements
  10.  
  11. * better documentation
  12.  
  13. Thomas Barrett
  14. Physics Dept, Ohio State University
  15. barrett@pacific.mps.ohio-state.edu
  16.  
  17. 18 Aug 93
  18.  
  19. // ------------------------- cut here -------------------------
  20.  
  21. This code is a hand-assembled version of the fft routine from Numerical
  22. Recipes.  See the book for information about how it works.  All variable
  23. names in comments refer to those in the book.
  24.  
  25. To use this routine:
  26.  
  27. * You must have a math coprocessor.
  28.  
  29. * Use Think C (users of other compilers may be able to adapt it).
  30.  
  31. * Set "Native floating-point format" under "Compiler Settings" in the
  32.     Options... dialog box.  This uses the 12-byte format which the math
  33.     coprocessor uses internally.
  34.  
  35. void tb_four1(long double *data, long nn, long isign);
  36.  
  37. * Store your data to be processed as an array of 12-byte long doubles.
  38.     Note that this will take more memory than the 8-byte doubles.  Also,
  39.     the array must be of 'nn' complex numbers, where each complex number
  40.     is a pair of long doubles.  'data' should therefore be a pointer to
  41.     an array of 24*nn bytes.
  42.  
  43. * This routine is DESTRUCTIVE.  The output data is placed in the space
  44.     where the input data was.  If you still want the input, make a copy
  45.     and pass the copy to the routine.
  46.  
  47. * In the book Numerical recipes in C, from which this routine is taken,
  48.     the first element of the array is accessed as data[1].  This is an
  49.     error!  C uses data[0] for the first element of an array.  In C, this
  50.     can be corrected by using data[i-1] and data[i] instead of data[i]
  51.     and data[i+1] (they always occur in pairs).  This routine expects
  52.     'data' to be a pointer to the first element of the array.  If you
  53.     are replacing the C version, and compensated for this in the routine
  54.     that called four1 (like the book suggest), then this is an issue.
  55.  
  56. * 'nn' must be a power of 2 (like 8, 16, 32...).  Useable range is between
  57.     8 and 128M (2^27 complex numbers).
  58.  
  59. * 'isign' must be 1 or -1, where -1 corresponds to an inverse fft.
  60.  
  61. * See the book for input and output formats.
  62.  
  63. I strongly recommend that, if you have an fft routine already working, you
  64. test this to make sure it gives the correct values when placed in your
  65. program (always a good idea).  It's been used successfully for a
  66. couple of months, and it is nearly twice as fast as the C version compiled
  67. by Think C 5 with optimizations.
  68.  
  69. Thomas Barrett
  70. Physics Dept, Ohio State University
  71. barrett@pacific.mps.ohio-state.edu
  72.  
  73. Thanks to:  Dan Flatin & Pascal Laubin.
  74.  
  75. #define    PI    3.1415926535897932384626433
  76.  
  77. /* ----------------------- tb_four1.c -----------------------------
  78.  
  79. optimized version of Numerical Recipes' fft routine
  80.  
  81. Thomas Barrett, 1993
  82.  
  83. written for Think C
  84. this routine assumes that data contains 12-byte 6888x-native long doubles
  85. also, you must have a math coprocessor to run this routine
  86.  
  87. -------------------------------------------------------------------
  88. register usage:
  89.  
  90. d0    I                a0    data[I]            fp0    WPR
  91. d1    J                a1    data[J]            fp1    WPI
  92. d2    M                a2    data[0]            fp2    WR
  93. d3    loop, MMAX                            fp3    WI
  94. d4    ISTEP                                fp4    TEMPR    \
  95. d5    NN,N                                fp5    TEMPI     \    internal
  96. d6    internal                            fp6             /    calculations
  97.                                         fp7            /
  98. ---------------------------------------------------------------- */
  99.  
  100. void tb_four1(long double *data, long nn, long isign)
  101. {
  102.     long double twopi = 2.0 * PI * isign;
  103.     
  104.     asm 68020, 68881 {
  105.         movem.l    a2/d3-d6,-(sp)
  106.         fmovem.x    fp4-fp7,-(sp)
  107.         
  108.         move.l    nn,d5
  109.         clr.l    d3
  110.         move.l    d5,d3        ; d3 = loop counter
  111.         move.l    #-1,d0        ; i(d0) = -1
  112.         movea.l    data,a0
  113.         suba.l    #12,a0        ; pointer to array indexed by 0
  114.         movea.l    a0,a2        ; a2 = *(data[0])
  115.         suba.l    #12,a0
  116.         
  117.     ; ------------ re-order values ---------------------------------
  118.         
  119.         move.l    #1,d1        ; j(d1) = 1
  120. @bits    adda.l    #24,a0        ; a0 = *(data[i])
  121.         addq.l    #2,d0        ; i += 2
  122.         cmp.l    d1,d0        ; cmp j,i
  123.         bge        @nosw        ; branch if i(d0) >= j(d1)
  124. @swap    movea.l    a2,a1
  125.         move.l    d1,d6
  126.     ;    mulu.l    #12,d6
  127.         lsl.l    #2,d6        ; these four instructions are equivalent to
  128.         adda.l    d6,a1        ; the mulu.l #12 and save a dozen cycles
  129.         lsl.l    #1,d6
  130.         adda.l    d6,a1        ; a1 = *(data[j])
  131.  
  132.         fmove.x    (a0),fp0    ; swap
  133.         fmove.x    (a1),fp1
  134.         fmove.x    fp1,(a0)
  135.         fmove.x    fp0,(a1)
  136.         fmove.x    12(a0),fp0
  137.         fmove.x    12(a1),fp1
  138.         fmove.x    fp1,12(a0)
  139.         fmove.x    fp0,12(a1)
  140.  
  141. @nosw    move.l    d5,d2        ; m(d2) = nn(d5) = #points
  142. @jloop    cmp.l    #2,d2
  143.         blt        @jrdy        ; branch if m(d2) < 2
  144.         cmp.l    d2,d1
  145.         ble        @jrdy        ; branch if j(d1) <= m(d2)
  146. @fixj    sub.l    d2,d1        ; j -= m
  147.         lsr.l    #1,d2        ; m /= 2
  148.         bra        @jloop
  149. @jrdy    add.l    d2,d1        ; j += m
  150.         subq.l    #1,d3
  151.         bne        @bits
  152.         
  153.     ; --------------- order is now ready -------------------------
  154.         
  155.         lsl.l    #1,d5        ; n(d5) = 2*nn(was d5) = #long doubles
  156.         move.l    #2,d3        ; mmax(d3) = 2
  157.         
  158.     ; -------------------- outer loop -----------------------------
  159.         
  160. @loop    cmp.l    d3,d5
  161.         ble        @done        ; branch if n(d5) <= mmax(d3)
  162.         move.l    d3,d4
  163.         lsl.l    #1,d4        ; istep(d4) = 2*mmax(d3)
  164.         fmove.x    twopi,fp1
  165.         fmove.l    d3,fp0
  166.         fdiv.x    fp0,fp1        ; theta(fp1) = 2 pi / mmax(d3)
  167.         fmove.x    fp1,fp0
  168.         fmove.w    #2,fp2
  169.         fdiv.x    fp2,fp0        ; fp0 = 1/2 theta
  170.         fsin.x    fp0
  171.         fmul.x    fp0,fp0
  172.         fmul.x    fp2,fp0
  173.         fneg.x    fp0            ; wpr(fp0) = -2 sin^2(1/2 theta)
  174.         fsin.x    fp1            ; wpi(fp1) = sin(theta)
  175.         fmove.w    #1,fp2        ; wr(fp2) = 1
  176.         fmove.w    #0,fp3        ; wi(fp3) = 0
  177.         
  178.     ; ------------------ inner loops -------------------------
  179.         
  180.         move.l    #1,d2        ; m(d2) = 1
  181. @mloop    move.l    d2,d0        ; i(d0) = m(d2)
  182.  
  183.         move.l    d0,d6        ; i(d0)
  184.         movea.l    a2,a0
  185.     ;    mulu.l    #12,d6
  186.         lsl.l    #2,d6
  187.         adda.l    d6,a0
  188.         lsl.l    #1,d6
  189.         adda.l    d6,a0        ; a0 = pointer to 1st i
  190.         movea.l    a0,a1
  191.         move.l    d3,d6        ; mmax(d3)
  192.     ;    mulu.l    #12,d6
  193.         lsl.l    #2,d6
  194.         adda.l    d6,a1
  195.         lsl.l    #1,d6
  196.         adda.l    d6,a1        ; a1 = pointer to 1st j
  197.         move.l    d4,d6        ; istep(d4)
  198.         mulu.l    #12,d6        ; 12 * istep. pointer increment
  199.  
  200. @iloop    move.l    d0,d1
  201.         add.l    d3,d1        ; j(d1) = i(d0) + mmax(d3)
  202.         
  203.     ;    movea.l    a2,a1
  204.     ;    move.l    d1,d6
  205.     ;    mulu.l    #12,d6
  206.     ;    adda.l    d6,a1        ; a1 = *(data[j(d1)])
  207.     ;    movea.l    a2,a0
  208.     ;    move.l    d0,d6
  209.     ;    mulu.l    #12,d6
  210.     ;    adda.l    d6,a0        ; a0 = *(data[i(d0)])
  211.  
  212.         fmove.x    (a1),fp4    ; fp4 = data[j]
  213.         fmove.x    fp4,fp7
  214.         fmul.x    fp2,fp4
  215.         fmove.x    12(a1),fp6    ; fp6 = data[j+1]
  216.         fmove.x    fp6,fp5
  217.         fmul.x    fp3,fp6
  218.         fsub.x    fp6,fp4        ; tempr(fp4) = wr(fp2)*data[j] - wi(fp3)*data[j+1]
  219.         fmul.x    fp2,fp5
  220.         fmul.x    fp3,fp7
  221.         fadd.x    fp7,fp5        ; tempi(fp5) = wr*data[j+1] + wi*data[j]
  222.         
  223.         fmove.x    (a0),fp6    ; fp6 = data[i]
  224.         fmove.x    fp6,fp7
  225.         fadd.x    fp4,fp6
  226.         fmove.x    fp6,(a0)    ; data[i] = data[i] + tempr(fp4)
  227.         fsub.x    fp4,fp7
  228.         fmove.x    fp7,(a1)    ; data[j] = data[i] - tempr(fp4)
  229.  
  230.         fmove.x    12(a0),fp6    ; fp6 = data[i+1]
  231.         fmove.x    fp6,fp7
  232.         fadd.x    fp5,fp6
  233.         fmove.x    fp6,12(a0)    ; data[i+1] = data[i+1] + tempi(fp5)
  234.         fsub.x    fp5,fp7
  235.         fmove.x    fp7,12(a1)    ; data[j+1] = data[i+1] - tempi(fp5)
  236.         
  237.         adda.l    d6,a0
  238.         adda.l    d6,a1
  239.  
  240.         add.l    d4,d0        ; i(d0) += istep(d4)
  241.         cmp.l    d5,d0
  242.         ble        @iloop        ; branch if i(d0) <= n(d5)
  243.         
  244.     ; ---------------update wr & wi ------------------------
  245.         
  246.         fmove.x    fp2,fp5        ; wtemp(fp5) = wr(fp2)
  247.         fmove.x    fp2,fp6
  248.         fmul.x    fp0,fp6
  249.         fadd.x    fp6,fp2        ; wr(fp2) += wr(fp2) * wpr(fp0)
  250.         fmove.x    fp3,fp6
  251.         fmul.x    fp1,fp6
  252.         fsub.x    fp6,fp2        ; wr(fp2) -= wi(fp3) * wpi(fp1)
  253.         fmove.x    fp3,fp6
  254.         fmul.x    fp0,fp6
  255.         fadd.x    fp6,fp3        ; wi(fp3) += wi(fp3) * wpr(fp0)
  256.         fmul.x    fp1,fp5
  257.         fadd.x    fp5,fp3        ; wi(fp3) += wtemp(fp5) * wpi(fp1)
  258.         
  259.         addq.l    #2,d2        ; m(d2) += 2
  260.         cmp.l    d3,d2
  261.         blt        @mloop        ; branch if m(d2) < mmax(d3)
  262.         move.l    d4,d3        ; mmax(d3) = istep(d4)
  263.         bra        @loop
  264.  
  265.     ; -------------------- done ----------------------------
  266.  
  267. @done    fmovem.x    (sp)+,fp4-fp7
  268.         movem.l    (sp)+,a2/d3-d6
  269.     }
  270. }
  271.